package com.siyu.model.xajmodel;


import cn.hutool.core.util.NumberUtil;

import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

public class XajModel {
    public static Map<Object,Object> calc(XajData xajData){
        /**以下为新安江输入参数获取**/
        //时间步长与步长类型
        int D=xajData.getPeriodTypeNum();
        //降雨及蒸发
        List<Double> pr = xajData.getPr();
        List<Double> er = xajData.getEr();
        //流域蓄水容量
        double WUM=xajData.getParamMap().get("wum");
        double WLM=xajData.getParamMap().get("wlm");
        double WDM=xajData.getParamMap().get("wdm");
        // 流域蒸散发折减系数
        double K = xajData.getParamMap().get("k");
        // 不透水面积比例(不透水面积占全面积百分数)
        double IMP =xajData.getParamMap().get("imp");
        // 流域蓄水容量不均匀分配系数
        double B = xajData.getParamMap().get("b");
        // 深层蒸散发系数
        double C = xajData.getParamMap().get("c");
        // 地下水消退系数
        double KKG = xajData.getParamMap().get("kkg");
        // 自由水蓄水容量（流域平均自由水蓄水容量，主要取决于下垫面的地质结构）
        double SM = xajData.getParamMap().get("sm");
        // 自由水容量曲线指数(反映流域自由水容量分布的均匀程度)
        double EX = xajData.getParamMap().get("ex");
        // 地下径流出流系数
        double KG = xajData.getParamMap().get("kg");
        // 壤中流出流系数
        double KSS = xajData.getParamMap().get("kss");
        // 壤中流消退系数
        double KKSS = xajData.getParamMap().get("kkss");

        /** 以下为新安江模型计算初始参数 **/
        // 上层土壤蓄水量
        double WU = 0;
        // 下层土壤蓄水量
        double WL = 70;
        // 深层土壤蓄水量
        double WD = 80;
        // 自由水蓄水量（初始流域自由水水深）
        double S = 20;
        // 初始产流面积比
        double FR = 0.1;
        // 壤中流初始值
        double QRI0 = 0.1;
        // 地下径流初始值
//        double QRG0 = 0.05;
        double QRG0 =xajData.getParamMap().get("jl");

        /**以下为计算及过程变量**/
        // 总径流流量
        List<Double> QR = new ArrayList<>();
        // 地面径流流量
        List<Double> QRS = new ArrayList<>();
        // 壤中流径流流量
        List<Double> QRI = new ArrayList<>();
        // 地下径流流量
        List<Double> QRG = new ArrayList<>();
        double rr=0;
        double ar=0;
        double pp=0;
        // 每等分的小时数
        // 流域土壤张力水容量
        double WM = WUM + WLM + WDM;
        // 步长内的壤中流蓄水库消退系数
        double CI = Math.pow(KKSS, 1.0 / D);// 公式换算
        // 步长内的地下水蓄水库消退系数
        double CG = Math.pow(KKG, 1.0 / D);// 公式换算
        // 步长内壤中流出流系数
        double KSSD = (1 - Math.pow(Math.abs(1 - (KSS + KG)), 1.0 / D)) / (1 + KG / KSS);// 公式换算
        // 步长内地下水出流系数
        double KGD = KSSD * KG / KSS;// 公式换算
        // 以时间长度进行循环计算
        for (int i = 0; i < pr.size(); i++) {
            double tempE = er.get(i);
            double tempP = pr.get(i);
            // 蒸发数据不能小于零
            if (tempE < 0)
                tempE = 0;
            // 降雨数据不能小于零
            if (tempP < 0)
                tempP = 0;
            // 扣除蒸发降雨量（主要变量）
            double PE = tempP - K * tempE;
            // 单元初始蓄水量计算，预热期前数据不准确
            double W0 = WU + WL + WD;
            // 蓄满产流量（主要变量）
            double R = 0.0;
            // 不透水面积产流量
            double RIMP = 0.0;
            if (PE >0) {
                // 张力水蓄水容量对应纵坐标
                double A = 0.0;
                // 流域最大点蓄水容量,透水面积上的最大蓄水容量
                double WMM = (1 + B) * WM / (1 - IMP);
                double WW = 1 - W0 / WM;
                if (WW < 0) {
                    WW = 0;
                }
                // 公式计算蓄水容量曲线横坐标
                A = (WMM * (1 - Math.pow(WW, 1 / (B + 1))));
                // if判断流域是否蓄满
                if (PE + A < WMM) {
                    // 未蓄满公式计算径流深
                    R = PE - WM + W0 + WM * Math.pow((1 - (PE + A) / WMM), 1 + B);
                } else {
                    R = PE - WM + W0;// 蓄满公式计算径流深
                }
                if (R < 0.0) {
                    R = 0.0;
                }
                if (R > PE) {
                    R = PE;
                }
                RIMP = PE * IMP;// 计算不透水面积上的产流量
            } else {
                double RG = S * KG * FR;
                double RI = RG * KSS / KG;
/*                if (R < 0) {
                    R = 0.0;
                }*/
                R=0;
                if (FR > 0) {
                    S = S - (RG + RI) / FR;//修改R为S
                }
            }
            // 更新上、下、深层张力水容量
            WU = WU + PE - R;
            // 当上层张力水容量小于零时
            if (WU < 0) {
                // 判断下层张力水容量
                if (WL <= C * WLM) {
                    WL = WL + C * WU;
                    WU = 0;
                    if (WL < 0) {
                        WD = WD + WL;
                        WL = 0;
                    }
                } else {
                    WL = WL + WU * WL / WLM;
                    WU = 0;
                    if (WL < 0) {
                        WD = WD + C * WL;
                        WL = 0;
                    }
                }
            }
            // 当上层张力水容量大于零量，需要判断是否大于最大蓄水容量
            else {
                // 当计算蓄水容量大于最大蓄水容量
                if (WU > WUM) {
                    // 当上层蓄水容量满后，则往下层传播
                    WL = WL + WU - WUM;
                    WU = WUM;
                    // 当下层蓄水容量大于最大蓄水容量
                    if (WL > WLM) {
                        // 当下层蓄水容量满后，则往深层传播
                        WD = WD + WL - WLM;
                        WL = WLM;
                    }
                }
            }
            if (WL < 0)
                WL = 0;
            if (WU < 0)
                WU = 0;
            if (WD < 0)
                WD = 0;
            // 三水源划分：
            // (1)重新划分步长
            // (2)各水源径流深计算

            // 5MM 地面径流
            double RS = 0.0;
            // 5mm 壤中流
            double RI = 0;
            // 5mm 地下径流
            double RG = 0;
            // 5mm 壤中流出流系数
            double KSSDD = 0;
            // 5mm 地下径流出流系数
            double KGDD = 0;

            // 三水源划分：
            // (1)重新划分步长
            // (2)各水源径流深计算
            // 5MM量级地面径流深
            if (PE <= 0) {
                RS = 0.0;
                RI = S * KSS * FR;
                RG = S * KG * FR;
                S = S - (RI + RG) / FR;
            } else {
                // 全流域最大点的自由水蓄水容量
                double SMM = (1 + EX) * SM; // 流域最大自由水容量
                double FR0 = FR;
                // 产流面积比
                FR = R / PE;
                // 自由水蓄水量计算
                S = S * FR0 / FR;
                double Q = R / FR;
                double NN;
                if (Q % 5 == 0.0) {
                    NN = Q / 5;
                } else {
                    NN = (int) (Q / 5) + 1;
                }
                Q = Q / NN;
                // 以5mm为尺度分割净雨
                KSSDD = ((1 - Math.pow((1 - (KSSD + KGD)), 1 / (double) NN)) / (1 + KGD / KSSD));// 公式转换
                KGDD = KSSDD * KGD / KSSD;// 公式转换
                // 产流面积上最大一点的自由水蓄水容量
                double SMMF = 0.0;
                if (EX <= 0.0001) // EX约为0，表明流域自由水容量分布均匀
                {
                    SMMF = SMM; // 产流面积上最大点自由水容量
                } else {
                    SMMF = (SMM * (1 - Math.pow(1 - FR, 1 / EX)));// 公式计算
                }
                // 产流面积上的自由水平均蓄水容量深
                double SMF = SMMF / (1 + EX); // 流域最大平均蓄水容量

                // 自由水蓄水容量对应坐标
                double AU = (SMMF * (1 - Math.pow((1 - S / SMF), 1 / (1 + EX))));// 公式计算
                for (int j = 0; j < NN; j++) {
                    if (S > SMF) {
                        S = SMF;
                    }
                    // 5MM 地面径流
                    double NRS = 0.0;
                    // 5mm 壤中流
                    double NRI = 0;
                    // 5mm 地下径流
                    double NRG = 0;
                    // 径流深公式计算
                    if (Q + AU <= 0) {
                        NRS = 0;
                        NRI = 0;
                        NRG = 0;
                        S = 0;
                    } else if (Q + AU >= SMMF) {// 判断是否需满
                        // 蓄满情况下各个水源净雨及自由水蓄水量计算
                        NRS = (Q + S - SMF) * FR;// 地面径流公式计算
                        NRI = SMF * KSSD * FR;// 壤中流径流公式计算
                        NRG = SMF * KGD * FR;// 地下径流公式计算
                        S = SMF - (NRI + NRG) / FR;// 自由水蓄水量公式计算
                    } else if (Q + AU < SMMF) {// 判断是否蓄满
                        // 未蓄满情况下的各个水源净雨计算
                        NRS = (Q - SMF + S + SMF * Math.pow(1 - (Q + AU) / SMMF, EX + 1)) * FR;// 地面径流公式计算
                        NRI = (S + Q - NRS / FR) * KSSDD * FR;// 壤中流流公式计算
                        NRG = (S + Q - NRS / FR) * KGDD * FR;// 地下径流公式计算
                        S = S + Q - (NRS + NRI + NRG) / FR;// 自由水蓄水量公式计算
                    }
                    RS += NRS;
                    RI += NRI;
                    RG += NRG;
                }
            }
            // 模型汇流计算：
            // (1)地面径流计算
            // (2)壤中流径流计算
            // (3)地下径流计算
            // 坡地汇流
            if (QRI0 < 0) {
                QRI0 = 0;
            }
            if (QRG0 < 0) {
                QRG0 = 0;
            }
            // 地面径流
//            QRS[i] = RS * (1 - IMP) + RIMP * IMP;
            QRS.add(RS  + RIMP);
            // 壤中流
            QRI.add(QRI0 * CI + RI * (1 - CI));
            // 地下径流
            QRG.add(QRG0 * CG + RG * (1 - CG));
            // 总出流
//            QR.add(QRS.get(i) + QRI.get(i)+QRG.get(i));
            QR.add(NumberUtil.round(QRS.get(i) + QRI.get(i)+QRG.get(i), 2).doubleValue());
            QRI0 = QRI.get(i);
            QRG0 = QRG.get(i);
            rr+=QR.get(i);
            pp+=tempP;

            if (FR < 0) {
                FR = 0.01;
            }
            if (S < 0) {
                S = 0;
            }
            R = 0;
            RG = 0;
            RS = 0;
        }
        if(pp == 0){
            ar=0.65;
        }else {
            ar=rr/pp;
        }
        if(ar>1.0){
            ar=0.65;
        }

        // 调用汇流计算
        double l = xajData.getParamMap().get("l");
        Double cs = xajData.getParamMap().get("cs");
        List<Double> SubOut = new ArrayList<>();// 出口流量序列
        List<Double> QW = new ArrayList<>();//流量过程
        double area = xajData.getArea();// 流域面积
        area=area/(xajData.getPeriodTypeNum()*3.6);//单位换算
        for(int i=0;i<QR.size();i++){
//            QW.add(QR.get(i)*area);
            QW.add(NumberUtil.round(QR.get(i)*area, 2).doubleValue());
        }
        for (int i =0; i < QW.size(); i++) {
            if(i<l||i<1){
                SubOut.add(0.0) ;
            }else{
                SubOut.add(NumberUtil.round(cs * SubOut.get(i - 1) + (1 - cs) * QW.get(i - (int)l), 2).doubleValue());
            }
        }
        HashMap<Object, Object> map = new HashMap<>();
        // 新安江计算的入库流量
        map.put("preQ",SubOut);
        // 径流系数
        map.put("runoffNum",ar);
        // 净雨量
        map.put("netRain",QR);
        return map;
    }
}
